Most data has a spatial dimension to it. Knowing ‘where’ the data is coming from is often as crucial as knowing the ‘what’, ‘when’ and ‘who’ dimensions of a given dataset. Therefore it should be no surprise that R has a rich suite of packages for constructing maps and analyzing spatial data. R’s capability has grown so much over the years that it’s functionality rivals many dedicated geographic information systems (GIS). During this Meetup the basics for managing and mapping spatial data with be introduced, using the following packages: sf, ggplot2, tmap, mapview and leaflet.
library(maps)
library(maptools)
map("state")
map("state", "indiana")
st <- map("state", fill = TRUE, plot = FALSE)
st_sp <- map2SpatialPolygons(st, IDs = st$names)
st_sp$state <- st$names
proj4string(st_sp) <- CRS("+init=epsg:4326")
# data frame
str(st)
## List of 4
## $ x : num [1:15599] -87.5 -87.5 -87.5 -87.5 -87.6 ...
## $ y : num [1:15599] 30.4 30.4 30.4 30.3 30.3 ...
## $ range: num [1:4] -124.7 -67 25.1 49.4
## $ names: chr [1:63] "alabama" "arizona" "arkansas" "california" ...
## - attr(*, "class")= chr "map"
# S3 object
str(st_sp, 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 63 obs. of 1 variable:
## ..@ polygons :List of 63
## ..@ plotOrder : int [1:63] 50 28 4 33 30 2 5 63 44 25 ...
## ..@ bbox : num [1:2, 1:2] -124.7 25.1 -67 49.4
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
plot(st_sp)
library(sf)
library(USAboundaries)
cnty <- us_counties()
cnty <- subset(cnty, !state_name %in% c("Alaska", "Hawaii"))
vars <- c("statefp", "countyfp")
plot(cnty[vars])
IN <- read_sf(dsn = "D:/geodata/soils/soils_GSMCLIP_mbr_2599033_03/wss_gsmsoil_IN_[2006-07-06]/spatial/gsmsoilmu_a_in.shp", layer = "gsmsoilmu_a_in")
# simplify polygons
IN <- rmapshaper::ms_simplify(IN)
library(aqp)
library(soilDB)
miami <- fetchKSSL("Miami")
idx <- complete.cases(miami$x, miami$y)
miami_sp <- miami[idx, ]
coordinates(miami_sp) <- ~ x + y
proj4string(miami_sp) <- CRS("+init=epsg:4326")
miami_sp <- SpatialPointsDataFrame(
coords = coordinates(miami_sp@sp),
data = site(miami_sp),
proj4string = CRS(proj4string(miami_sp))
)
miami_sf <- st_as_sf(miami_sp)
miami_sf <- within(miami_sf, {
lon = st_coordinates(miami_sf)[, 1]
lat = st_coordinates(miami_sf)[, 2]
})
https://github.com/potterzot/rnassqs
# devtools::install_github('potterzot/rnassqs')
# library(rnassqs)
# source("C:/Users/steph/Nextcloud/code/api_keys.R")
data(state)
st <- state.abb
# corn_us <- lapply(st, function(x) {
# cat("getting", x, as.character(Sys.time()), "\n")
# tryCatch({
# corn = nassqs_yield(
# list("commodity_desc"="CORN",
# "agg_level_desc"="COUNTY",
# "state_alpha"=x
# ),
# key = nass_key
# )},
# error = function(err) {
# print(paste("Error occured: ",err))
# return(NULL)
# }
# )
# })
# corn_us <- do.call("rbind", corn_us)
#
# save(corn_us, file = "C:/workspace2/corn_us.RData")
# write.csv(corn_us, file = "nass_corn_us.csv", row.names = FALSE)
load(file = "C:/Users/Stephen.Roecker/Nextcloud/data/corn_us.RData")
corn_yield <- subset(corn_us, short_desc == "CORN, GRAIN - YIELD, MEASURED IN BU / ACRE")
corn_yield <- within(corn_yield, {
Value = as.numeric(Value)
year = as.numeric(year)
state_name = NULL
})
cnty_corn <- merge(cnty, corn_yield,
by.x = c("state_abbr", "countyfp"),
by.y = c("state_alpha", "county_code"),
all.x = TRUE
)
corn_states <- c("IL", "IA", "IN", "MI", "MN", "MO", "NE", "OH", "SD", "ND", "WI")
library(dplyr)
library(ggplot2)
group_by(corn_yield, state_alpha, year) %>%
summarize(
yield_low = min(Value, na.rm = TRUE),
yield_median = median(Value, na.rm = TRUE),
yield_max = max(Value, na.rm = TRUE)
) %>%
filter(state_alpha %in% corn_states) %>%
mutate(state = state_alpha,
source = "NASS"
) %>%
ggplot() +
geom_line(aes(x = year, y = yield_median, col = source)) +
geom_ribbon(aes(x = year, ymin = yield_low, ymax = yield_max), alpha = 0.25) +
facet_wrap(~ state) +
ylab("yield per county (bu/acre)") +
ggtitle("USDA-NASS Corn Yields")
# geom_point(data = yld_sum[yld_sum$state %in% corn_states, ], aes(x = 2018, y = yield_med, col = "NASIS"), size = 1) +
# geom_ribbon(data = yld_sum2[yld_sum2$state %in% corn_states, ], aes(x = year, ymin = yield_low2, ymax = yield_max2, col = "NASIS"), alpha = 0.25) +
# geom_pointrange(data = yld_sum[yld_sum$state %in% corn_states, ], aes(x = 2018, y = yield_med, ymin = yield_low2, ymax = yield_max2, col = source))
library(ggplot2)
# loads data from maps package
st <- map_data("state")
# tidy example data
st_tidy <- broom::tidy(st_sp, region = "state")
IN_tidy <- broom::tidy(as(IN, "Spatial"), region = "MUSYM")
dim(st_sp)
## [1] 63 1
dim(st)
## [1] 15537 6
# Lines
ggplot() +
geom_point(data = miami_sf, aes(x = lon, y = lat)) +
geom_path(data = st, aes(x = long, y = lat, group = group)) +
xlim(range(miami_sf$lon)) +
ylim(range(miami_sf$lat)) +
ggtitle("Location of Miami Lab Pedons")
# Polygons
ggplot() +
geom_polygon(data = st_tidy, aes(x = long, y = lat, group = group, fill = id)) +
coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
# remove legend
guides(fill = FALSE)
# sf package
# Polygons
ggplot() +
geom_sf(data = cnty, aes(fill = statefp, lty = NA)) +
geom_sf(data = miami_sf) +
coord_sf(crs = "+init=epsg:5070") +
guides(fill = FALSE)
# Facets
test <- cnty_corn[cnty_corn$year %in% 2012:2017, ]
ggplot() +
geom_sf(data = test, aes(fill = Value, lty = NA)) +
scale_fill_viridis_c(na.value = "transparent") +
facet_wrap(~ year) +
geom_path(data = st, aes(x = long, y = lat, group = group)) +
ggtitle(corn_yield$short_desc[1])
library(ggmap)
# get ggmap
bb <- sf::st_bbox(IN)
bb <- make_bbox(lon = bb[c(3, 1)], lat = bb[c(2, 4)])
gmap <- get_map(bb, maptype = "terrain", source = "osm")
# Lines
ggmap(gmap) +
geom_path(data = IN_tidy, aes(x = long, y = lat, group = group))
# geom_sf(data = IN, fill = NA, inherit.aes = FALSE) +
# guides(fill = FALSE)
# geom_sf() doesn't work with ggmp, their is a systematic shift https://github.com/r-spatial/sf/issues/336
library(tmap)
tm_shape(IN) + tm_polygons("MUSYM", border.col = NULL) +
tm_shape(cnty) + tm_borders() +
tm_shape(miami_sf) + tm_dots() +
tm_legend(legend.outside = TRUE)
# interactive web map
tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(IN) + tm_borders()
library(mapview)
cols <- RColorBrewer::brewer.pal(50, "Paired")
test <- mapview(IN, zcol = "MUSYM", lwd = 0, col.regions = cols) +
mapview(cnty, type = "l")
test
# export to html
mapshot(test, url = "C:/workspace2/test.html", selfcontained = FALSE)
library(leaflet)
test <- leaflet() %>%
addProviderTiles("Esri.WorldImagery", group = "Imagery") %>%
addPolygons(data = IN, fill = FALSE, color = "black", weight = 2)
test
# export to html
htmlwidgets::saveWidget(test, file = "C:/workspace2/test.html", selfcontained = FALSE)
Healy, K., 2018. Data Visualization: a practical introduction. Princeton University Press. http://socviz.co/
Gimond, M., 2019. Intro to GIS and Spatial Analysis. https://mgimond.github.io/Spatial/
Hijmans, R.J., 2019. Spatial Data Science with R. https://rspatial.org/
Lovelace, R., J. Nowosad, and J. Muenchow, 2019. Geocomputation with R. CRC Press. https://bookdown.org/robinlovelace/geocompr/
Pebesma, E., and R. Bivand, 2019. Spatial Data Science. https://keen-swartz-3146c4.netlify.com/